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ABSTRACT 


In  order  to  increase  building  safety  under  earthquake  motions,  there  has  been 
increasing  interest  in  base  isolation  with  passive  isolators.  Computer  modeling  is  an 
important  aspect  of  the  building  design  and  evaluation  process,  but  solving  for  the 
transient  response  of  large  structural  systems  with  localized  nonlinearities  is 
computationally  demanding.  Current  finite  element  programs  can  rapidly  detennine 
normalized  mode  shapes  and  natural  frequencies  of  several  thousand  degree  of  freedom 
structures  for  use  in  determining  the  transient  response.  However,  actual  computation  of 
the  transient  response  can  be  very  time-consuming  and  expensive  for  such  large 
structures.  A  recently  developed  convolution  algorithm  utilizes  the  Volterra  integral  in  a 
recursive  block-by-block  integral  equation  fonnulation  to  efficiently  compute  the 
transient  response  of  multi-story,  nonlinear,  base  isolated  buildings.  This  algorithm  was 
utilized  in  a  versatile  optimization  scheme  which  determines  parameters  for  both  linear 
and  nonlinear  mathematical  model  isolators  coupled  to  a  multi-degree  of  freedom 
structure.  To  optimize  the  isolator  parameters,  the  procedure  incorporates  modal 
properties  computed  from  a  finite  element  model  of  the  structure,  the  earthquake 
accelogram  of  interest,  and  user-defined  objective  and  constraint  functions.  An  example 
is  given  of  a  4-story,  single  bay  structure  subjected  to  the  1940  El  Centro  NS  excitation. 
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I.  INTRODUCTION 


The  use  of  earthquake  protective  systems  has  increased  greatly  in  the  last  30  years 
with  numerous  structures  around  the  world  being  built  or  retrofitted  with  seismic 
isolation  systems.  Seismic  isolation,  or  base  isolation,  is  a  design  strategy  based  on  the 
premise  that  it  is  feasible  to  uncouple  a  structure  from  the  ground  in  order  to  protect  it 
from  the  damaging  effects  of  earthquake  motions  [Ref.  1].  The  base  isolation  system 
provides  additional  flexibility  and  damping  designed  to  absorb  the  earthquake  energy  and 
thereby  reduce  the  severity  of  earthquake  attacks  on  the  structure.  A  basic  feature  of  base 
isolation  systems  is  that  the  system  restricts  large  defonnations  to  special  components, 
the  isolators,  while  the  building  vibrates  almost  as  a  rigid  body  [Ref.  2]. 

The  design  process  for  a  structural  base  isolation  system  is  very  challenging  due 
to  the  inherent  randomness  of  earthquakes  as  well  as  regional  characteristics  such  as 
geology,  soil  composition,  proximity  to  faults,  and  potential  wind  loading  effects.  The 
primary  goal,  of  course,  is  to  design  the  system  to  best  protect  the  structure  and  its 
contents  from  earthquake  attack.  In  this  regard,  numerical  optimization  is  an  effective 
tool  in  obtaining  suitable  parameters  for  use  in  preliminary  design  efforts,  but  dynamic 
analysis  of  large,  complex  structures  is  also  computationally  demanding.  This  large 
computational  cost  strongly  inhibits  the  iterative  procedures  required  in  traditional 
optimization  methods. 

However,  more  efficient  techniques  have  been  developed  to  calculate  structural 
responses.  A  recently  developed  convolution  algorithm  utilizes  modal  infonnation  from 
a  finite  element  model  and  the  Volterra  integral  in  a  recursive  block-by-block  integral 
equation  formulation  (RBBIF)  to  efficiently  compute  the  transient  response  of  multi¬ 
story,  nonlinear,  base  isolated  buildings  [Ref.  3].  This  algorithm  is  extremely  fast  and 
was  utilized  in  a  versatile  optimization  scheme  capable  of  detennining  parameters  for  a 
library  of  linear  and  nonlinear  mathematical  model  isolators  coupled  to  a  multi-degree  of 
freedom  (DOF)  structure.  Within  the  optimization  procedure,  the  RBBIF  algorithm 
makes  use  of  time  and  frequency  domain  synthesis  techniques  to  rapidly  recalculate 
structural  system  responses  of  the  finite  element  model  as  the  parameter  design  variables 
are  altered  with  each  iteration.  To  optimize  the  isolator  parameters,  the  algorithmic 
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procedure  incorporates  modal  properties  computed  from  the  finite  element  model  of  the 
structure,  the  earthquake  accelogram  of  interest,  and  user-defined  objective  and  constraint 
functions. 


A.  ALGORITHM  OVERVIEW 

The  optimization  algorithm  described  in  this  paper  is  very  versatile  and  quite 
suitable  for  use  in  dynamic  analysis  of  structures  with  localized  nonlinearities  [Ref.  3]. 
Base  isolation  systems  essentially  uncouple  the  structure  from  the  ground,  so  by 
exploiting  the  preprocessing  capabilities  of  the  program,  any  unconstrained,  free-free 
finite  element  model  can  be  used.  To  determine  the  transient  response  using  RBBIF, 
modal  properties  are  first  extracted  from  the  free-free  finite  element  structure  for  use  in 
later  function  calls. 

The  program  can  detennine  numerous  items  of  interest  to  designers  and  planners 
such  as  base  displacement,  top  floor  acceleration,  relative  velocity,  maximum  values  for 
displacement,  velocity,  and  acceleration,  and  root  mean  square  values  for  displacement, 
velocity,  and  acceleration.  In  addition,  the  algorithm  can  be  easily  modified  to  determine 
other  items  of  interest. 

There  is  also  a  library  of  linear  and  nonlinear  mathematical  isolator  models.  Of 
frequent  use  in  current  passive  isolation  design  are  the  bilinear  and  Wen  hysteretic 
models.  The  uplift  isolator  is  modeled  as  a  linear  spring  per  [Ref.  4],  but,  nonlinearities 
can  also  be  incorporated  into  this  portion  of  the  program.  The  library  of  isolators  is 
discussed  further  in  the  paper  with  an  emphasis  on  the  bilinear  and  Wen  elements. 

It  is  possible  to  input  any  design  earthquake  relevant  to  a  particular  region  of 
interest.  Since  the  program  is  very  fast,  the  potential  exists  to  input  different  design 
earthquakes  to  detennine  if  similar  isolator  parameter  values  are  obtained.  Earthquake 
excitations  can  also  be  scaled  within  the  program.  Numerous  earthquake  time  histories 
are  available  at  websites  such  as  [Ref.  5]  and  [Ref.  6].  The  design  earthquake  used  in  the 
later  example  is  the  El  Centro  1940  North-South  (NS)  component. 
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II.  EARTHQUAKES 


Events  such  as  volcanic  activity,  explosions,  and  collapsing  cave  roofs  may  cause 
seismic  activity.  However,  the  most  important  earthquakes  from  an  engineering 
standpoint  are  of  tectonic  origin,  or,  in  other  words,  those  associated  with  large-scale 
strains  in  the  crust  of  the  earth  [Ref.  7].  Numerous  analyses  associate  slip  along  geologic 
faults  as  a  primary  mechanism  to  produce  tectonic  motions,  particularly  in  areas  laced 
with  many  faults  such  as  Southern  California.  Recent  tectonic  earthquakes  have  occurred 
with  devastating  consequences  in  Turkey,  India,  Iraq,  Japan,  and  California. 

Earthquakes  of  magnitude  5.0  or  greater  on  the  Richter  scale  generate  ground 
motions  severe  enough  to  cause  significant  structural  damage  [Ref.  8].  Regarding 
structural  design,  prediction  of  the  applied  seismic  intensity  may  be  complicated  and 
filled  with  uncertainties.  Structures  on  competent  bedrock  are  subjected  predominantly 
to  the  effects  of  relatively  high-frequency,  low-amplitude  vibrations  for  relatively  short 
durations.  The  ground  surface  in  such  settings  is  not  likely  to  suffer  pennanent 
deformation  or  ground  failure  during  an  earthquake.  Structures  on  compressible  deposits, 
however,  particularly  where  the  water  table  is  high  (within  33  feet  of  the  surface),  are 
subjected  not  only  to  the  effects  of  relatively  low  frequency,  high  amplitude  vibrations, 
but  also  may  be  subjected  to  disruption  caused  by  differential  settlement,  lateral 
displacement,  or  liquefaction  [Ref.  9].  This  was  seen  in  the  damaged  structures  of  the 
Marina  District  in  San  Francisco  after  the  1989  Loma  Prieta  earthquake. 


A.  EARTHQUAKE  CHARACTERISTICS 

Because  earthquake  motions  are  irregular  and  each  earthquake  is  different  from 
all  others,  even  at  a  given  site,  it  is  important  to  establish  whatever  characteristics  certain 
groups  of  earthquakes  may  have  in  common  and  to  then  base  earthquake -resistant  design 
on  these  generalizations  [Ref.  7].  [Ref.  7]  describes  four  general  groups  of  earthquakes: 

1)  Practically  a  single  shock.  Figure  1  is  the  East-West  (EW)  component  of  the 
earthquake  which  occurred  in  Port  Hueneme,  California,  on  March  18,  1957. 
These  occur  only  at  short  distances  from  the  epicenter,  only  on  firm  ground,  and 
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only  for  shallow  earthquakes.  Single  shock  earthquakes  are  characterized  with 
moderate  magnitudes  (5. 4-6.2),  shallow  foci  (less  than  19  miles),  almost 
unidirectional  motion,  and  a  prevalence  of  short  periods  of  vibration  (0.2  seconds 
or  less).  This  type  of  earthquake  motion  is  very  uncommon. 
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Figure  1.  Port  Hueneme  1957  EW  Component  (From:  [Ref.  7,  p.226]) 


2)  Moderately  long,  extremely  irregular  motion.  Figure  2  is  the  North-South  (NS) 
component  of  the  El  Centro  earthquake  which  occurred  on  May  18,  1940.  These 
earthquakes  occur  at  a  moderate  distance  from  the  epicenter  and  only  on  firm 
ground.  Moderately  long  and  irregular  earthquakes  are  characterized  by  a  wide 
range  of  vibration  periods  (0.05-0.5  sec,  and  2.5-6  sec)  and,  ordinarily,  equal 
severity  in  all  directions  rather  than  the  unidirectional  motion  associated  with 
single  shock  earthquakes.  Of  note,  almost  all  earthquakes  originating  along  the 
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Circumpacific  Belt  are  of  this  type,  so  design  earthquakes  such  as  the  El  Centro 
1940  NS  component  can  be  applicable  for  similar  geological  regions. 
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Figure  2.  El  Centro  1940  NS  Component  (From:  [Ref.  7,  p.  227]) 


3)  Long  ground  motion  with  pronounced  prevailing  periods  of  vibration  (e.g.,  any 
Mexico  City  earthquake).  Figure  3  is  the  NS  component  of  the  earthquake  which 
occurred  in  Mexico  City  on  July  6,  1964.  This  type  of  earthquake  motion  results 
from  filtering  either  (1)  a  single  shock  earthquake,  or  (2)  a  moderately  long 
earthquake  with  irregular  motion,  through  layers  of  soft  soil  within  the  range  of 
linear  or  almost  linear  soil  behavior.  These  earthquake  motions  are  further 
amplified  by  successive  wave  reflections  at  the  interfaces  of  these  soil  mantles. 
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Figure  3.  Mexico  City  1964  NS  Component  (From:  [Ref.  7,  p.  227]) 


4)  Ground  motion  involving  very  large-scale,  pennanent  defonnations  of  the  ground, 
with  possible  slides  or  soil  liquefaction  (e.g.,  Anchorage  earthquake  of  1964). 
Ordinarily,  it  is  impractical  to  attempt  a  structural  design  to  resist  large-scale 
failure  of  the  ground.  In  regions  where  seismic  studies  indicate  this  type  of 
motion  is  probable,  it  is  best  to  erect  the  structure  at  a  different  site  or  somehow 
treat  the  soil  in  such  a  way  that  the  phenomenon  becomes  unlikely,  at  least 
locally. 

There  are,  of  course,  ground  motions  with  characteristics  spanning  each  of  these 
broad  categories,  but  this  classification  serves  the  purpose  of  outlining  general 
characteristics  which  may  be  useful  in  base  isolation  design. 


B.  DESIGN  EARTHQUAKES 

General  seismic  factors  which  must  be  considered  are  local  geology  and  soil 
properties,  proximity  to  faults,  recorded  histories  of  earthquakes  in  the  region,  and  known 
or  probable  characteristics  of  earthquakes  in  the  region  such  as  intensity  or  period  [Ref. 

6 


10].  The  program  works  well  in  these  respects  since  any  scaled  or  actual  earthquake  time 
history  can  be  applied  to  the  structural  finite  element  model. 

The  design  earthquakes  available  from  [Refs.  5  and  6]  are  downloaded  in  the  form 
of  data  points.  Typical  time  steps  are  on  the  order  of  10'  ,  so  each  earthquake  time 
history  has  hundreds  of  data  points.  The  MATLAB  interp  function  is  used  to  interpolate 
between  the  points  in  order  to  obtain  a  complete  time  history  for  use  in  the  program. 

With  knowledge  of  the  seismic  factors  in  the  region  of  interest,  applicable  design 
earthquake  time  histories  from  the  previously  mentioned  website  resources  can  be 
obtained  and  applied  in  the  optimization  program  to  detennine  whether  the  isolator 
design  parameters  remain  consistent  and  suitable  for  each  regional  earthquake. 


C.  HORIZONTAL  AND  VERTICAL  RESPONSES 

The  destructive  motion  in  an  earthquake  is  usually  the  horizontal  ground 
movement  [Ref.  11].  Although  large  vertical  components  have  been  recorded  in  some 
recent  earthquakes,  it  is  felt  that  these  will  not  usually,  of  themselves,  be  destructive  but 
are  only  a  problem  if  coupling  occurs  between  them  and  the  horizontal  components. 
Such  coupling  is  avoided  in  the  system  described  in  [Ref.  1 1]  by  manufacturing  bearings 
with  large  vertical  stiffness  and  low  horizontal  stiffness,  taking  care  to  ensure  the 
fundamental  period  of  the  base  isolated  building  is  greater  than  the  periods  of  the  regional 
design  earthquakes. 

Although  the  majority  of  the  destructive  ground  motion  tends  to  be  horizontal,  the 
RBBIF  algorithm  can  efficiently  compute  both  horizontal  and  vertical  structural 
responses.  The  development  of  the  structural  synthesis  is  discussed  in  subsequent 
sections,  and,  in  the  later  example,  a  4-story  finite  element  structure  is  subjected,  in  the 
horizontal  direction,  to  the  frequently  used  El  Centro  1940  NS  component  design 
earthquake. 
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III.  LIBRARY  OF  SEISMIC  ISOLATOR  MODELS 


As  previously  mentioned,  a  library  of  isolator  models  is  available  for  use  in  the 
program.  The  library  consists  of  linear  and  nonlinear  mathematical  models  representing: 

1 .  Linear  spring 

2.  Linear  spring  and  viscous  damper  in  parallel 

3.  Friction  isolator 

4.  Ideal  bilinear  element  with  constant  yield  stress  properties 

5.  Real  bilinear  element  with  changing  yield  stress  properties 

6.  Maxwell  element 

7.  Wen  element 

Of  frequent  use  in  base  isolation  design  are  the  bilinear  element  and  the  Wen 
element.  These  represent  frequency-independent  isolators  which  are  capable  of 
maintaining  a  unique  hysteretic  loop  across  a  wide  range  of  frequencies.  The  Maxwell 
element  is  also  frequency-independent,  but  it  is  not  studied  in  this  thesis. 

Material  properties  can  vary,  but,  passive  isolators  known  as  laminated-rubber 
bearings  can  generally  be  modeled  in  a  linear  manner  since  these  have  a  linear  restoring 
force  and  linear  damping  [Ref.  12].  Typical  nonlinear  passive  isolators  which  could  be 
modeled  by  the  bilinear,  Wen,  or  Maxwell  elements  are  high-damping  rubber  bearings 
(HDRB),  elastomeric  lead-rubber  bearings,  friction  dampers,  steel  dampers,  and  lead- 
extrusion  dampers. 


A.  HYSTERETIC  ISOLATORS 

The  bilinear  element  and  the  Wen  element  are  part  of  a  group  collectively  referred 
to  as  hysteretic  isolators.  The  term,  hysteretic,  is  regarded  as  frequency-independent 
damping  associated  with  Coulomb  friction  at  the  material  grain  boundaries,  and  it  also 
refers  to  the  offset  between  the  loading  and  unloading  curves  under  cyclic  loading  or 
earthquake  excitation  [Ref.  7].  The  bilinear  and  Wen  hysteretic  models  have  gained 
recognition  as  accurate  and  useful  tools  to  numerically  portray  various  structural  damping 
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behaviors.  The  general  nonlinear  equation  that  describes  frequency  independent  isolators 
is: 


F=-K0(l  +  iS)(x-ug)  (1) 

where: 


K0  _  dynamic  stiffness 
8  -  described  as  the  loss  factor 
x  -  base  motion 
ug  -  base  motion 

Hysteretic  energy  dissipation  is  one  of  the  most  effective  means  of  providing  the 
substantial  levels  of  damping  required  of  a  base  isolation  system.  Generally  speaking, 
large  displacements  can  be  effectively  controlled  if  additional  substantial  damping  is 
introduced  into  the  structure  through  the  base  isolators.  Figure  4  shows  an  idealized 
bilinear  force-displacement  loop  where  the  enclosed  area  is  a  measure  of  the  energy 
dissipated  during  one  cycle  of  motion.  Ki  and  Kj  are  two  distinct  stiffness  constants,  and 
yp  is  the  yield-to-post  ratio  which  relates  Kj  and  K2.  The  bilinear  model  is  an 
approximation  to  the  hysteretic  loop  shown  in  Figure  5. 
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Tensile 


Figure  4.  Ideal  Bilinear  Force-Displacement  Loop 


In  general  tenns,  Equation  (1)  describes  the  hysteretic  loop  shown  in  Figure  5. 
The  mechanism  which  forms  the  hysteretic  loop  is  cyclical  deformation  along  a 
material’s  stress-strain  curve.  Deformation  below  the  tensile  and  compressive  yield 
stresses  describes  typical  linear  elastic  behavior.  To  begin  generating  the  hysteresis  loop, 
deformation  occurs  beyond  the  initial  yield  stress,  resulting  in  plastic  deformation.  In  the 
plastic  region,  the  material  resists  additional  stress  through  strain-hardening,  which  leads 
to  a  secondary  stiffness  rate  (AT  in  Figure  4).  Unloading  from  the  plastic  region  occurs 
elastically  until  plastic  deformation  occurs  again  as  the  cyclic  loading  increases. 
Unloading  from  the  plastic  region  then  occurs  elastically  again,  which  completes  the 
hysteretic  loop  and  begins  the  next  cycle. 
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Restoring 

force 


Figure  5.  Hysteretic  Loop 


Examples  of  rubber  isolators  which  utilize  hysteretic  energy  dissipation  to  absorb 
earthquake  excitations  are  the  high-damping  rubber  bearings  and  elastomeric  lead-rubber 
bearings  mentioned  previously  [Ref.  12]. 

Rubber  bearings  such  as  these  are  widely  used  in  passive  isolation  designs.  In 
their  basic  form,  rubber  bearings  provide  vertical  support,  horizontal  flexibility,  and 
centering  forces.  By  inserting  additional  materials  such  as  carbon  black  or  lead  plugs,  the 
hysteretic  energy  dissipation  feature  is  increased  as  damping  is  increased  [Refs.  1 1  and 
12]. 
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B.  WEN  ELEMENT  MODEL 

The  Wen  element  is  a  mathematical  equation  which  accurately  models  the  hysteretic 
behavior  of  the  previously  mentioned  nonlinear  rubber  isolators.  The  ensuing  description  is 
summarized  from  [Ref.  2]. 

The  Wen  nonlinear  equation  refines  the  bilinear  hysteretic  loop  from  Figure  4  to 
generate  a  more  realistic  hysteresis  by  varying  the  parameters  of  the  isolator  restoring  force. 
Equation  (2): 


m  =  ccf  x(t)  +  (1  -  a)Fyz(t) 


(2) 


where: 


OC—x(t)  —  equivalent  linear  portion 
Sr 

(1  —  ci)F  z(t')=  the  nonlinear  portion 

CC  —  post-to-preyielding  stiffness  ratio 

8y  —  yield  displacement  of  isolator 

Fy  =  yield  force  of  isolator 

z(t)  =  dimensionless  hysteretic  displacement 


z(t)  is  defined  by  the  nonlinear  first  order  differential  equation: 


Sv  z  =  -y 


ZZ^-Pxz'  +  Ay 


(3) 


where: 


(3,  y,  A  =  dimensionless  parameters  which  control  the  shape  of  the  hysteresis  curve 


//  =  integer  which  controls  smoothness  of  transition  from  elastic  to  plastic 
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The  combination  of  the  nonlinear  parameters  from  Equation  (3)  yields  the  maximum 
value  of  the  restoring  force. 
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IV.  RECURSIVE  BLOCK-BY-BLOCK  FORMULATION 


The  recently  developed  RBBIF  convolution  algorithm  rapidly  computes  structural 
transient  responses  of  a  finite  element  model  [Ref.  3].  This  rapid  analysis  capability 
facilitated  development  of  an  efficient  and  versatile  numerical  optimization  scheme  for 
the  design  of  nonlinear,  hysteretic  isolators.  An  overview  of  the  RBBIF  formulation  is 
presented  here.  The  complete  fonnulation  is  covered  in  [Refs.  3,  14,  and  15]. 


A.  BACKGROUND 

Generating  the  equations  of  motion  of  an  N-degree  of  freedom  (NDOF)  structure 
allows  the  structure  to  be  described  mathematically  both  in  a  static  and  a  dynamic 
manner.  The  number  of  DOF  is  directly  controlled  when  modeling  the  structure  in  a 
finite  element  program.  By  controlling  the  number  of  DOF,  the  number  of  differential 
equations  describing  the  structure  is  then  controlled  and,  hence,  computationally 
manageable.  In  order  to  discretize  a  structure  with  an  infinite  number  of  DOF,  a  logical 
determination  must  first  be  made  with  regards  to  retaining  a  sufficient  number  of  DOF 
which  sufficiently  describe  the  structure. 

Following  this  step,  the  structure  can  then  be  described  by  a  continuous 
mathematical  model  and  formulated  as  a  set  of  partial  differential  equations  (PDE).  The 
PDEs  can  then  be  transformed  into  a  set  of  2nd  order  ordinary  differential  equation  (ODE) 
by  using  elements  of  lumped-parameter  models.  The  fonnulation  is  covered  in  structural 
dynamics  textbooks  such  as  [Ref.  8]. 

Structural  synthesis  refers  to  substructure  coupling  and  structural  modification  in 
the  finite  element  model.  In  the  context  of  the  physical  coordinate  synthesis  formulations 
developed  in  [Refs.  14  and  15],  a  structural  system  is  defined  to  consist  of  one  or  more 
uncoupled  substructures.  A  single  governing  equation  for  nonlinear  transient  synthesis 
was  derived  in  [Ref.  3]  which  addressed  each  of  the  following  three  general  analysis 
category  sets: 

1.  Structural  Modification  (m-set)  -  the  addition  and/or  removal  of  linear  and/or 
nonlinear  structural  elements 
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2.  Prescribed  Base  Motion  (b-set)  -  application  of  base  motion  to  structure 
through  linear  and/or  nonlinear  elements 

3.  Substructure  Coupling  (c-set)  -  the  joining  of  substructures  (a  linear  analysis) 

The  “m-sef “b-set”,  and  “c-set”  are  subsets  of  DOF  associated  with  structural 
modification  (such  as  base  isolators),  prescribed  base  motion,  and  substructure  coupling, 
respectively.  Each  subset  may  include  DOF  from  all  substructures.  As  reported  in  [Ref. 
14],  the  nonlinear  elements  are  installed  in  the  synthesis.  Consequently,  all  the 
substructures  are  linear,  and  coupling  is  also  a  linear  synthesis. 

Substructure  coupling  allows  nonlinear  elements  to  be  isolated  by  division  of  the 
system  into  substructures.  This  feature  of  the  program  results  in  significant 
computational  savings  since  the  eigensolutions  of  the  linear  substructures  are  only  solved 
once.  Within  the  optimization  program,  the  coupling  set,  or  c-set,  is  of  primary  interest 
since  this  defines  the  set  of  DOF  subjected  to  the  nonlinear  earthquake  excitation. 

Each  substructure  is  described  by  impulse  response  functions  calculated  at  the 
DOF  where  nonlinear  elements  are  to  be  installed,  where  loads  are  applied,  and  at  other 
DOF  for  which  synthesized  nonlinear  transient  response  is  required. 

For  the  linear  substructures,  the  IRF  are  most  efficiently  calculated  using  modal 
superposition  as  described  in  [Ref.  8].  In  order  to  decrease  computational  cost,  however, 
the  RBBIF  algorithm  only  utilizes  a  sufficient  number  of  modes  to  ensure  convergence  of 
the  IRF  within  a  specified  tolerance.  These  converged  IRF  are  virtually  indistinguishable 
from  the  “exact”  IRF,  as  shown  in  [Ref.  14]. 


B.  GOVERNING  EQUATIONS 

The  discretized  finite  element  structural  equation  of  motion,  assuming 
proportional  structural  damping,  is: 

+tct~{}i  +  [KL„UL,  ={CL  (5) 

where: 

n  =  total  number  of  DOF  in  the  system 
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{  Fc°eJ(t)  }  =  external  isolator  restoring  forces 

The  general  solution  to  this  second  order  ordinary  differential  equation  is  the  sum 
of  the  homogeneous  and  particular  solutions: 

x(t)  =  x(t)h  +  x{t)p  (6) 

The  homogeneous  solution  for  a  specific  DOF  is: 

x(t)h  =  e~C6>"'(A  sin  yJl~C B  cos  V1  -  C MJ)  (7) 

The  derivation  of  the  particular,  or  transient,  solution  depends  upon  the  form  of 
the  forcing  function.  For  linear  forces,  such  as  those  of  a  linear  spring,  the  solution  is 
attainable  through  manipulation  of  Equation  (5).  However,  for  nonlinear  elements  such 
as  hysteretic  isolators,  the  particular  solution  can  be  very  complex.  The  particular 
solution  can  be  obtained  by  using  the  governing  equation  for  transient  structural 
synthesis,  a  nonlinear  Volterra  integral  equation  involving  a  convolution-type  kernel. 
The  derivation  of  this  convolution  integral  is  developed  in  mechanical  vibration 
textbooks  such  as  [Ref.  16],  and  is  formulated  as: 

x(t)p=\‘_h(t-T)f(T)dT  (8) 

where: 

f(t)  =  design  earthquake  time  history 

h(t)  =  matrix  of  impulse  response  functions  (IRF) 


C.  TIME  DOMAIN  SYNTHESIS 

Substituting  the  particular  solution,  Equation  (8),  into  Equation  (6)  results  in  the 
total  solution  for  the  general  equation  of  motion  written  in  expanded  format  here: 


Xx(t) 

xpo 

x2(t ) 

: 

— 

X2(t) 

s1  ... 

Xm(I) 

+ 


W-T) 

hjt-z) 


olA.fi  “D 


hJt-T) 

hJt-T) 


KAt-r) 


where: 


KH-r) 

FJt-T) 


f;\t) 


■dr 


(9) 


hJt-T) 


0 
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m  =  number  of  DOF  in  the  coupling  set 
n  =  number  of  modes 


The  {  F’so{t )  }  vector  components  represents  the  isolator  restoring  forces  resulting  from 

ground  excitation  on  the  nodes  denoted  by  the  subscript. 

To  reduce  the  computational  requirements  for  Equation  (9),  Equation  (5)  is 
transformed  into  modal  space.  The  linear  modal  transformation  equation  is  defined  as: 


where: 

<f>/  =  mode  shapes 
n  =  number  of  modes 
i  =  ith  mode 
j  =jth  node 


j= i 


(10) 


By  substituting  Equation  (10)  into  Equation  (5),  and  then  pre-multip lying  with  the 
transpose  of  the  mass  nonnalized  mode  shape  matrix,  [<f>]T  ,  the  modal  space  equivalent 
becomes: 


<hO) 

2^°l 

0 

0 

. 

?i(0 

0  • 

•  0 

?1« 

[*! 

*2  ' 

Fx(t) 

+ 

0 

2^2w2 

0 

?2(0 

+ 

0 

CO2  ■ 

•  0 

q2(t) 

= 

•i  ■ 

'  •*] 

F2(t) 

”„(0 

0 

0 

2^ncon 

q„  (0 

0 

0  • 

• 

*2  • 

■ 

Fn«\ 

(11) 


The  modal  solution  of  Equation  (11)  yields  the  equivalent  Volterra  integral  form 
of  Equation  (12): 


?!<<) 

9l(0 

rt 

q2(0 

= 

?m(0 

h 

0 

hn(t-T)  0 

0  h22(t-T) 


0 

0 

hnn(t-T) 


[4>r 


1J 

•2] 

4>::1 


\dr 


(12) 
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Since  the  system  of  equations  was  transformed  into  modal  space  the  impulse 
response  matrix,  [h(t-x)],  is  uncoupled  and  denoted  by  the  (~)  symbol.  Equation  (12)  is 
then  transformed  back  into  nonnal  space  by  multiplying  with  [(f)]  and  using  the  modal 
transformation  relationship,  Equation  (10).  This  results  in: 


where: 

m  =  number  of  DOF  in  the  coupling  subset 
n  =  number  of  modes 

F‘so  =  DOF  in  the  coupling  subset  which  are  subjected  to  the  excitation  force. 

Within  the  RBBIF  algorithm,  the  E]wo(f)  components  correspond  to  the  seismic  isolators. 

From  Equation  (13),  it  is  seen  that  there  are  significantly  more  unknowns  than 
equations.  There  are  (n+m)  unknowns  and  only  (n)  coupled  equations.  However,  since 
there  are  (m)  forces  in  {/7,'"(f)} ,  a  recursive  iteration  process  can  be  used  on  a  reduced 
system  of  coupled  equations  which  is  (m  x  m)  in  size.  This  total  solution  is  written  as: 


As  previously  mentioned,  the  RBBIF  algorithm  only  retains  a  specified  number  of 
modes,  based  on  a  user-defined  cutoff  frequency,  which  further  reduces  the 
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computational  cost.  By  utilizing  this  feature  and  only  solving  for  the  DOF  in  excitation 
(c-set),  the  overall  result  is  a  large  savings  in  time  and  computational  effort. 


D.  RECURSIVE  ITERATION 

To  further  decrease  the  computational  time  and  effort,  a  recursive  iteration 
formulation  is  used.  To  illustrate  the  recursive  iteration  process,  an  example  is  given  of  a 
nonlinear  structure,  Figure  6,  subjected  to  excitation,  y(t).  The  nonlinear  structure  is 
isolated  by  two  linear  springs  with  respective  stiffness  coefficients,  Kj  and  K2. 


Figure  6.  Nonlinear  Structure  Isolated  by  Two  Linear  Springs  with  Kj,  K2 


The  isolator  restoring  forces  are: 

fi(t)  =Kj{  x,(t)-y(t)  } 

(15) 

fi(t)  =  K2{  x2(t)-y(t)  } 

(16) 

Or  in  matrix  fonn: 

Kx  0  U(0-y(0j 

0  K2\x2(t)-y(t) } 


{/(')}  = 


(17) 
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By  using  Equation  (8),  the  Volterra  convolution  integral,  to  solve  for  the 
particular  solution,  the  transient  response  appears  as: 


V  o  jvo -y(t)\dT 

0  K2  \x2(t)-y(t)\ 

where  H(t-x)  is  known: 


{x(t)}  =  \tH(t-T) 


(18) 


H(t-r) 


Hn(t-T )  0 

H22(t-T) 


which  results  in  Equation  (19)  after  substituting  Equation  (18): 


VO 


VO 


ff 0  ipv 

H21(t-T)\  l/2(0j 


(19) 


(20) 


Since  {x(()}  and  {f(t)  |  are  unknown  vectors,  arbitrary  values  for  components  of 
{ x(t)}  are  chosen,  and  the  solution  is  determined  with  an  iterative  process  until  a  user- 
specified  convergence  tolerance  is  met.  The  iteration  procedure  is  as  follows: 

1.  Select  {x(t)}.  Set  {x0id(t)}  =  {x(t)} . 

2.  Use  {x0id(t)}  to  compute  { fnew0 )}  from  Equation  (17). 

3.  Solve  for  {xnew(t)}  with  Equation  (20),  using  the  {fnew(t)}  calculated  in  the 
previous  step. 

4.  Check  if  {xnew(t)}  -  {x0id(t)}  <  convergence  tolerance. 

5.  If  convergence  tolerance  is  not  met,  set  {x0id(t)}  =  {xnew(t)}. 


6.  Repeat  steps  2-5  until  {xnew(t)}  -  {x0id(t)}  <  convergence  tolerance. 
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E.  RECURSIVE  BLOCK-BY-BLOCK  SYNTHESIS 

Incorporation  of  a  recursive  block-by-block  time  history  synthesis  further 
decreases  the  computational  time.  To  begin  with,  the  entire  time  history  is  divided  into 
blocks  of  time,  and  a  small  time  step  is  chosen  for  use  in  the  convolution  integral  to 
determine  the  time  history  response  over  the  period  of  each  block.  The  blocks  take 
advantage  of  the  linear  properties  of  the  convolution  integral,  and  the  time  history 
responses  computed  for  each  block  form  the  complete  response  over  the  entire  time 
history. 

As  a  general  example,  the  time  history  from  0  to  tn  is  divided  into  (n)  blocks,  0-ti, 
ti-t2,  to-t3,  ...,  tn.i-tn.  The  recursive  iteration  procedure  described  in  the  previous  section 
is  then  perfonned  for  each  block,  and  the  complete  solution  is  of  the  form: 

{'"  H(t-T){f(z)}dT  =  \l  H (7,  -r){/(T)}r/T  +  ['2  H(t2 -T){f{T)}dz+... 


-+r  H(tB-T){f(T)}dT 

Jtn- 1 

For  a  single  DOF,  Equation  (2 1)  is  written  in  matrix  form  as: 


m' 

x(A  t) 


'kll 
<  k}  > 
k). 


k(h)  J 

x(tt  +  At) 
-  : 
x(t2) 

x(t2  +  At) 
x(t3) 


[HA  [0]  [0] 

[HA  [H2]  [0] 

[HA  [HA  [HA 

[HA  [HA  [HA 


[0] 

[0] 

[0] 

[HA 


m 

/(AO 

m 

f(t\+  At) 

‘  f(t2)  > 

f(t2+  At) 

f(h) 


no 


where  Hi. . . //„  are  lower  triangular. 


(21) 


(22) 
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To  simplify,  let: 


fa  =  the  forcing  function  for  the  first  block,  where  fa  =  {/'(()) 

fb  =  the  forcing  function  for  the  second  block,  where  //,  =  \f(ti) .../(tf)} 


f„  =  the  forcing  function  for  the  nth  block,  where  f,  =  {f(tn-i)  —f(tn)} 


The  block-by-block  synthesis  is  completed  by  carrying  out  the  matrix 
multiplication  and  applying  the  recursive  iteration  in  the  following  manner: 

1.  Use  recursive  iteration  procedure  for  {xa}  and  \f,  }  until  convergence  is 
reached  for  {xa}  =  [Hi]  {fa} . 

2.  Use  recursive  iteration  procedure  for  {xb}  and  {//,]  until  convergence  is 
reached  for  Jx/,}  =  [If] \fa\  +  [//_?]{//,},  where  [If] {fa}  is  now  known  from 
step  (1). 

3.  Use  recursive  iteration  procedure  for  {xc}  and  {/',  }  until  convergence  is 
reached  for  {xc}  =  [If]  [f, \+[ff]  [fb\+[Hi]  iff,  where  [///]{/„}  is  known  from 
step  (1),  and  [/U]  {fb}  is  now  known  from  step  (2). 

4.  Continue  until  recursive  block-by-block  synthesis  is  complete. 


For  most  efficient  use  of  the  RBBIF  algorithm,  the  first  data  point  within  each 
block  must  be  correct.  It  is  necessary  that  the  force  solution  for  each  block  have  this 
correct  initial  value  in  order  to  obtain  an  accurate  time  history  between  blocks.  For  linear 
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isolators  such  as  the  two  springs  in  the  example,  the  recursive  block-by-block  iterations 
are  relatively  straightforward  between  the  individual  blocks. 

For  nonlinear  hysteretic  isolators,  however,  an  accurate  RBBIF  synthesis  relies 
heavily  on  force  solutions  from  the  previously  solved  block.  In  order  to  achieve 
continuity  of  the  hysteretic  isolator  response  between  time  blocks,  the  initial  block  value 
is  interpolated  with  a  direct  quadratic  approximation  from  the  last  three  force  solution 
data  points  of  the  previous  block  shown  in  Figure  7.  A  MATLAB  Runge-Kutta  4  (RK4), 
ordinary  differential  equation  solver  is  used  for  the  interpolation  of  Fpred  in  Figure  7. 
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Figure  7.  Recursive  Interpolation 
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V.  OPTIMIZATION 


Base  isolation  may  be  used  to  provide  effective  solutions  for  a  wide  range  of 
seismic  design  problems.  Buildings  such  as  hospitals,  fire  stations,  and  police  stations 
must  remain  in  operation  at  all  times,  particularly  when  a  disaster  such  as  a  major 
earthquake  has  caused  significant  devastation.  It  is  now  generally  accepted  that  base 
isolated  buildings  will  perform  better  than  conventional  fixed  base  buildings  in  moderate 
or  strong  earthquakes  [Ref.  17].  There  is  also  a  strong  economic  incentive  to  incorporate 
base  isolation  systems  into  building  design.  In  addition  to  minimizing  damage  from  an 
earthquake,  base  isolation  systems  may  also  lower  the  overall  cost  of  a  structure  in  terms 
of  less  material  requirements  in  the  design  and  decreased  insurance  rates.  Successful 
implementation  of  these  systems  could  potentially  result  in  enormous  savings  of  both 
lives  and  money,  but  the  design  process  is  challenging. 

Numerical  optimization  can  be  a  useful  tool  in  the  design  process  of  base 
isolators,  particularly  when  making  evaluations  during  the  preliminary  design  phase.  The 
RBBIF  method  of  determining  transient  structural  response  is  extremely  fast  and 
computationally  efficient,  so  its  use  in  an  optimization  scheme  is  very  appealing. 


A.  GENERAL  PROBLEM  FORMULATION 

In  an  optimization  problem,  the  quantity  to  be  maximized  or  minimized  is  the 
objective  function,  and  the  variables  to  be  altered  in  order  to  optimize  the  objective 
function  are  termed  design  variables.  The  constraint  functions  are  also  functions  of  the 
design  variables,  and  these  can  be  classified  as  inequality  constraints,  equality 
constraints,  or  side  constraints.  In  the  constrained  optimization  problem  solved  by  this 
program,  the  user  defines  the  objective  function  and  the  constraint  functions.  The 
objective  function  is  minimized,  while  ensuring  the  constraint  functions  are  not  violated 
within  a  user-defined  tolerance.  The  goal  in  the  formulation  of  the  optimization 
computer  programs  is  to  provide  a  means  of  determining  optimal  values  for  linear  and 
nonlinear  isolator  parameters,  in  order  to  minimize  a  dynamic  response  such  as  base 
displacement  or  top  floor  acceleration,  while  simultaneously  satisfying  all  prescribed 
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constraints.  Two  minimization  problems  were  studied  in  this  thesis.  The  first  was  to 
minimize  the  base  displacement,  and  the  second  was  to  minimize  the  top  floor 
acceleration.  The  general  problem  formulations  follow. 


1.  First  Minimization  Problem 

The  first  general  optimization  problem  formulation  in  the  program  is: 


Minimize:  (objective  function) 
/  ( parameters )  = 


[  ||**Be(0|Ll 

>-< 

max  xB  t(t) 

\<i<n 

l  Vnax  J 

Vnax 

(23) 


Ratio  of  Base  Displacement  °o-norm  to  Maximum  Structural  Displacement 


Subject  to:  (constraints) 


vR  (Oil  <0.70*X 

Base\  y ||oo  max 

Max  Base  Displacement  <  0.70*Initial  Max  Structural  Displacement 
Min  Limit  <  Isolator  Parameter  Values  <  Max  Limit 


(24) 


Design  Variables:  Isolator  Model  Parameters 


2,  Second  Minimization  Problem 

The  second  general  optimization  problem  formulation  in  the  program  is: 
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Minimize:  (objective  function) 


/  ( parameters )  = 


■^TopFloor^f)  ^ 

k  -  < 

max  xT  Fl  (t) 

1  <i<n 

Vnax 

'  -  s 

Vnax 

(25) 


Ratio  of  Top  Floor  Acceleration  co-norm  to  Maximum  Structural  Acceleration 


Subject  to:  (constraints) 

||Ww-(0|L<°.75*xmax 

(26) 

||^(0|L^0.70*ximx 

(27) 

Max  Top  Floor  Acceleration  <  0.75*Initial  Max  Structural  Acceleration 
Max  Base  Displacement  <  0 . 7 0 '"Initial  Max  Structural  Displacement 
Min  Limit  <  Isolator  Parameter  Values  <  Max  Limit 


Design  Variables:  Isolator  Model  Parameters 


B.  MAXIMUM  STRUCTURAL  DISPLACEMENT 

Determining  the  maximum  response  of  a  structure  to  all  possible  environmental 
excitations  is  a  fundamental  task  in  structural  dynamics,  but  the  maximum  response 
within  a  structure  may  vary  from  location  to  location  [Ref.  18].  A  method  reported  in 
[Ref.  19]  derives  operator  norms  for  the  convolution  operators  associated  with  linear, 
time-invariant  systems.  This  method  obtains  the  absolute  magnitude  of  the  Euclidean  2 
or  oo-norm  for  the  time  domain  response  of  a  multi-output  system  to  certain  classes  of 
input  disturbance  [Ref.  19].  This  method  was  further  refined  in  [Ref.  18]  to  determine 
the  maximum  response  of  a  structure  to  an  uncertain,  vector-valued  input,  such  as  an 
earthquake  time  history. 

Since  the  RBBIF  algorithm  is  capable  of  computing  nonlinear  dynamic  response 
through  use  of  the  Volterra  convolution  integral,  results  from  [Refs.  18  and  19]  were 
incorporated  in  order  to  efficiently  determine  the  maximum  absolute  response  of  the 
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entire  structure.  The  complete  formulation  is  found  in  [Refs.  18  and  19],  and  the  portion 
used  within  the  algorithm  is  summarized  here. 


1.  Linear  Operator  Norms 

The  convolution  integral  shown  in  Equation  (28)  is  a  linear  operator  acting  on 
{f(t)\  to  produce  {x(t)}\ 


x(t)  =  J  h(t  —  T)f \r)dr 


(28) 


where: 


x(t)  =  vector  of  response  outputs 


f(t)  =  applied  input  design  earthquake 


h(t)  =  matrix  of  impulse  responses  with  each  impulse  response  function  defined  as: 


hit) 


x(t) 


<  .  c. 


V  V 


x(  0)  +^g0)  gin  ^ cog 


(29) 


J) 


[Ref.  19]  provides  the  mathematical  basis  for  defining  the  norm  of  a  vector¬ 
valued  function  of  time.  For  an  n-dimensional  real-valued  vector  function  { x(t)},  [Ref. 
18]  defines  the  r-norm  to  be: 


|*(0|| 


1m 


i=l 


(30) 


Normally,  only  the  2-nonn  or  the  °°-nonn  are  of  interest.  When  r  is  2,  Equation  (30)  is 
the  Euclidean  norm.  When  r  is  this  is  the  maximum  absolute  value  component  of  the 
vector. 

By  extending  the  definition  of  a  vector  norm  to  consider  not  only  the  norm  across 
components  as  a  function  of  time,  but  the  norm  across  time  as  well,  [Ref.  18]  defines  the 
(p,  r)  norm  of  the  vector-valued  function  {/(/)  }  to  be: 
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(31) 


LDKolIJ  dt 


Of  note,  Equation  (31)  does  not  vary  with  time.  If  r  approaches  2  and  p 
approaches  2,  this  norm  is  the  root  mean  square  magnitude.  If  r  approaches  °°  and  p 
approaches  °°,  this  nonn  is  the  maximum  value,  over  all  time,  of  the  maximum  value  of 
any  component  in  the  vector  { x(t)}  [Ref.  18]: 

xmax=|HL00=  sup  |x(0|L  =  sup  max  1^(01  (32) 

’  —00<t<00  —oo<t<oo 


Equation  (31)  was  incorporated  into  the  optimization  algorithm  to  detennine  the 
maximum  structural  displacement,  velocity,  and  acceleration  over  all  time. 


C.  MATLAB  OPTIMIZATION  TOOLBOX 

The  optimization  tool  used  in  the  algorithm  is  the  MATLAB  fmincon  function, 
which  is  designed  for  optimization  of  a  constrained,  nonlinear,  multivariable  problem. 
Specific,  user-oriented  instructions  are  found  in  the  MATLAB  Optimization  Toolbox, 
[Ref.  20].  In  a  constrained  optimization  problem,  the  goal  is  to  transfonn  the  problem 
into  an  unconstrained  subproblem  that  can  then  be  solved  and  used  as  the  basis  of  an 
iterative  process.  The  basis  for  the  subproblem  is  satisfying  the  Kuhn-Tucker  (KT) 
conditions,  which  are  necessary  conditions  for  a  point,  x,  to  be  a  local  minimum.  The 
KT  conditions  are  stated  as: 

V/(x  ) -  ■  Vgj(x) -  •  Vhk  (x)  =  0  (33) 

j= 1  k=\ 

where: 

f(x)  =  objective  function 

g(x)  =  inequality  constraint  functions 

h(x)  =  equality  constraint  functions 

Uj  =  Lagrange  multipliers  >  0 

va  =  Lagrange  multipliers,  sign  unrestricted 
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The  Lagrange  multipliers,  u  and  v,  are  unspecified  parameters  which  convert  the 
constrained  problem  into  an  unconstrained  problem.  Only  active  constraints  are  included 
in  Equation  (33),  so  the  inactive  constraints  are  associated  with  Lagrange  multipliers 
equal  to  zero.  The  solution  of  Equation  (33)  forms  the  basis  to  many  nonlinear 
programming  algorithms  commonly  referred  to  as  Sequential  Quadratic  Programming 
(SQP)  methods.  The  finincon  function  uses  SQP  methods  to  detennine  a  search  direction 
for  the  design  variables  at  every  iteration.  An  overview  of  the  SQP  method  is  found  in 
[Ref.  21],  and  the  optimization  algorithm  flowchart  is  shown  in  Figure  8. 


Figure  8.  Optimization  Algorithm  Flowchart 
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VI.  FINITE  ELEMENT  BUILDING  MODEL  USED  IN  EXAMPLE 


In  order  to  demonstrate  the  capabilities  of  the  program,  a  finite  element  model  of 
a  4-story,  square  base,  single  bay  structure  was  used.  The  building  is  25  feet  wide  and 
has  an  inter-story  height  of  17.5  feet.  The  columns  and  beams  are  50  ksi  steel  structural 
members  with  the  following  specifications:  columns,  W36x486;  first  and  second  floor 
beams,  W36xl70;  third  floor  beams,  W36xl60;  fourth  floor  beams,  W36xl50;  roof 
beams,  W36xl35.  Each  floor  was  designed  for  a  150  psf  dead  load,  and  the  roof  was 
designed  for  a  50  psf  dead  load.  The  building  is  modeled  with  20  nodes,  120  DOF,  and 
four  DOF  in  excitation  in  the  horizontal  direction.  A  frame  illustration  of  the  building  is 
shown  in  Figure  9. 


Figure  9.  Single  Bay,  4-Story  Building  Frame 
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The  isolators  chosen  for  the  example  problems  are  real  bilinear  element  models 
with  the  following  initial  properties,  except  where  noted: 

Elastic  Stiffness:  k  =  15000  lbf/in 
Yield-to-Post  Ratio:  yp  =  0. 10 
Maximum  Tensile  Force:  fen  =  25000  lbf 
Maximum  Compressive  Force:  fcom  =  30000  lbf. 

The  building  weight  per  isolator  is  7740  lbf/isolator.  The  excitation  input  is  the  El  Centro 
1940  NS  earthquake  component,  and  a  cutoff  frequency  of  2  Hz  was  chosen  leading  to  12 
retained  modes  out  of  60  modes  total. 
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VII.  EXAMPLES  AND  RESULTS 


Four  examples  are  presented  in  which  optimal  bilinear  model  parameters  are 
recommended  for  the  4-story  structure  subjected,  in  the  horizontal  direction,  to  the  El 
Centro  1940  NS  component  design  earthquake.  A  Dell  Dimension  4100  computer  was 
used  to  run  the  optimization  program  with  a  time  step  of  0.05  sec.  The  primary  objective 
of  the  first  three  examples  was  to  reduce  the  maximum  base  displacement  to  70%  of  the 
maximum  structural  displacement,  xmax.  The  initial  stiffness  value,  k,  was  modified  in 
each  example.  The  primary  objective  of  the  last  example  was  to  reduce  the  maximum  top 
floor  acceleration  to  75%  of  the  maximum  structural  acceleration,  xmax .  The  initial 

parameters  for  the  bilinear  isolator  model  were: 

Examples  1  and  4:  k  =  15000  lbf/in 
yp  =  0.10 
fen  =  25000  lbf/in 

fcom  =  30000  lbf/in 

Example  2:  k  =  5000  lbf/in 

yp  =  0.10 
ften  =  25000  lbf/in 

fcom  =  30000  lbf/in 

Example  3:  k  =  20000  lbf/in 

yp  =  0.10 
fen  =  25000  lbf/in 

fcom  =  30000  lbf/in 

A.  GENERAL  PROBLEM  FORMULATIONS 

1.  First  Example:  Minimize  Base  Displacement,  k  =  15000  lbf/in 


33 


The  first  optimization  problem  was  to  minimize  the  peak  base  displacement  to 
less  than  70%  of  the  maximum  structural  displacement.  The  general  optimization 
formulation  was: 


Design  Variables:  Bilinear  Model  Parameters  ( k ,  y,„  ften,  fCOm) 


Minimize: 


/  ( parameters )  = 


f  !*&»  (OIL 

max  xBasei(t) 

\<i<n 

[  Vnax 

Vnax 

Subject  to:  g,  =  ^(OL  ^  0.70  *xmax 

g2  =  6000  <  k  <  25000 
g3  =0.01  <yp  <0.40 
g4  =  5000  <  ften  <  35000 
g5  =  -35000  <fcom<  -5000 


(34) 

(35) 

(36) 

(37) 

(38) 

(39) 


The  initial  parameters  for  the  bilinear  model  were: 
k=  15000  lbf/in 
yP  =  0.10 
fen  =  25000  lbf/in 
fcom  =  30000  lbf/in 


2.  Results  of  First  Minimization  Example 

The  optimization  program  recommended  the  following  bilinear  model  parameters 
in  order  to  minimize  the  base  displacement  response  to  less  than  70%  of  the  initial 
maximum  structural  displacement: 

k=  13463  lbf/in 
yP  =  0.273 
fen  =  25392  lbf 
fom  =  29544  lbf 
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The  minimized  objective  function  value  was  0.674,  or,  in  other  words,  the  base 
displacement  was  reduced  to  67.4%  of  the  initial  maximum  structural  displacement. 

The  optimization  program  accomplished  this  in  33  iterations.  Each  iteration  was 
approximately  136  seconds,  so  the  program  accomplished  this  within  75  minutes. 

Figure  10  shows  the  base  displacements  for  the  initial  and  the  optimized  design 
parameters.  The  oscillation  amplitude  is  increased  on  the  optimized  design  even  though 
the  overall  maximum  displacement  was  decreased.  The  increased  amplitude  of  the 
optimized  design  is  probably  because  the  yield-to-post  ratio,  yp,  increased,  which  made 
the  secondary  stiffness,  K2,  larger  than  the  initial  design  K2. 

Figures  11  and  12  show  hysteresis  loops  for  the  initial  and  optimized  designs, 
respectively.  The  optimized  hysteresis  loop  has  much  more  area  enclosed  than  the  initial 
design.  This  indicates  much  more  of  the  earthquake  energy  was  absorbed,  resulting  in 
decreased  base  displacement  magnitude. 


Initial  Design  -  Base  Isolated  Displacement 


Figure  10.  Base  Displacement  for  Initial  and  Optimized  Bilinear  Model  Parameters 
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Force  (Ibf)  «  Force  (Ibf) 


3.  Second  Example:  Minimize  Base  Displacement,  k  =  5000  Ibf/in 

The  second  optimization  problem  was  identical  to  the  first:  minimize  peak  base 
displacement  to  less  than  70%  of  the  maximum  structural  displacement. 

The  initial  parameters  for  the  bilinear  model  were: 
k  =  5000  lbf/in 
yP  =  o.io 

fen  =  25000  lbf/in 
fcom  =  30000  lbf/in 


4.  Results  of  Second  Minimization  Example 

The  optimization  program  recommended  the  following  bilinear  model  parameters 
in  order  to  minimize  the  base  displacement  response  to  less  than  70%  of  the  initial 
maximum  structural  displacement: 

k  =  6000  lbf/in 
yP  =  0.01 
fen  =  24740  lbf 
fcom  =  30004  lbf 

The  minimized  objective  function  value  was  0.646,  or  the  optimized  base 
displacement  was  reduced  to  67.4%  of  the  initial  maximum  structural  displacement. 

The  optimization  program  accomplished  this  in  only  1 1  iterations.  Each  iteration 
was  approximately  99  seconds,  so  the  program  accomplished  this  within  18  minutes. 

Figure  13  shows  the  base  displacements  for  the  initial  and  the  optimized  design 
parameters.  More  oscillations  are  apparent,  and  the  oscillation  amplitude  is  again 
increased  on  the  optimized  design  even  though  the  overall  maximum  displacement  was 
decreased.  This  is  probably  because  the  stiffness  value  as  well  as  the  yield-to-post  ratio 
are  small,  resulting  in  a  smaller,  more  flexible  secondary  stiffness,  AT 

Figures  14  and  15  show  hysteresis  loops  for  the  initial  and  optimized  designs, 
respectively. 
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Hysteresis  of  Base  Isolator  -  Optimized  Design 


Figure  15.  Hysteresis  Loop  for  Optimized  Bilinear  Model  Parameters 


5.  Third  Example:  Minimize  Base  Displacement,  k  =  20000  Ibf/in 

The  third  optimization  problem  was  also  identical  to  the  first:  minimize  peak 
base  displacement  to  less  than  70%  of  the  maximum  structural  displacement. 

The  initial  parameters  for  the  bilinear  model  were: 
k  =  20000  lbf/in 
yP  =  0.10 

ften  =  25000  lbf/in 
fcom  =  30000  lbf/in 


6.  Results  of  Third  Minimization  Example 

The  optimization  program  recommended  the  following  bilinear  model  parameters 
in  order  to  minimize  the  base  displacement  response  to  less  than  70%  of  the  initial 
maximum  structural  displacement: 

£  =  23271  lbf/in 
yP  =  0.24 
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fen  =  26434  lbf 
fcom  =  28733  lbf 


The  minimized  objective  function  value  was  0.697,  or  the  optimized  base 
displacement  was  reduced  to  69.7%  of  the  initial  maximum  structural  displacement. 

The  optimization  program  accomplished  this  in  39  iterations.  Each  iteration  was 
approximately  174  seconds,  so  the  program  accomplished  this  within  1 14  minutes. 

Figure  16  shows  the  base  displacements  for  the  initial  and  the  optimized  design 
parameters.  The  large  and  rapid  base  displacement  oscillations  seen  in  the  previous 
examples  are  not  apparent  here.  However,  the  optimized  design  top  floor  accelerations 
seen  in  Figure  17  are  larger  than  the  accelerations  from  the  initial  design. 

Figures  18  and  19  show  hysteresis  loops  for  the  initial  and  optimized  designs, 
respectively. 


Initial  Design  -  Base  Isolated  Displacement 


Figure  16.  Base  Displacement  for  Initial  and  Optimized  Bilinear  Model  Parameters 
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Figure  17.  Top  Floor  Acceleration  for  Initial  and  Optimized  Bilinear  Model  Parameters 
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Hysteresis  of  Base  Isolator  -  Initial  Design 


0  12  3 

Isolator  Displacement  (in) 


Figure  18.  Hysteresis  Loop  for  Initial  Bilinear  Model  Parameters 


Hysteresis  of  Base  Isolator  -  Optimized  Design 


-1  0  1 
Isolator  Displacement  (in) 


Figure  19.  Hysteresis  Loop  for  Optimized  Bilinear  Model  Parameters 
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Fourth  Example:  Minimize  Top  Floor  Acceleration 


The  second  optimization  problem  was  to  minimize  the  peak  top  floor  acceleration 
to  less  than  75%  of  the  maximum  structural  acceleration.  The  constraint  functions  from 
the  previous  example  were  maintained,  so  peak  base  displacement  remained  constrained 
to  70%  or  less  of  the  maximum  structural  displacement.  The  general  optimization 
formulation  was: 

Design  Variables:  Bilinear  Model  Parameters  (k,  yp,  fen,  fCOm) 


Minimize: 


/  ( parameters )  = 


TopFloor  (0  ^ 

>  -  < 

max  xT  Fl  (t) 

1  <i<n 

Vnax 

'  -  s 

Xmax 

(40) 


Subject  to:  S\  Ipro^FZoorCOl^  —  0-75  *xmax 

Si  =  I^SoseCOlL  -  °-70*Xmax 
g3  =  6000  <  k  <  25000 
g4  =0.01  <yp  <0.40 
g5  =  5000  <  fen  <  35000 
g6  =  -35000  <fcom<  -5000 


(41) 

(42) 

(43) 

(44) 

(45) 

(46) 


8.  Results  of  Fourth  Minimization  Problem 

The  optimization  program  recommended  the  following  bilinear  model  parameters 
in  order  to  minimize  the  top  floor  acceleration  response  to  less  than  75%  of  the  initial 
maximum  structural  acceleration: 

k=  13096  lbf/in 

yP  =  0.01 

fen  =  23232  lbf/in 
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fcom  =  22467  lbf/in 


The  minimized  objective  function  value  was  0.747.  The  top  floor  acceleration 
was  reduced  to  74.7%  of  the  initial  maximum  structural  acceleration. 

The  optimization  program  accomplished  this  in  36  iterations.  Each  iteration  was 
approximately  136  seconds,  so  the  program  accomplished  this  within  82  minutes. 

Figure  20  shows  the  top  floor  accelerations  for  the  initial  and  the  optimized 
design  parameters.  Although  the  maximum  top  floor  acceleration  magnitude  was 
decreased  by  the  optimized  design,  the  oscillation  amplitude  is  much  larger  over  the  latter 
portions  of  the  time  history.  This  indicates  that  larger  inertial  forces  would  occur  in  the 
structure  than  would  have  otherwise  occurred  with  the  initial  design,  even  though  the 
maximum  acceleration  is  minimized  by  the  optimized  design  parameters.  This  may  also 
indicate  the  optimized  design  parameters  are  not  really  good,  and  perhaps  a  different 
objective  function  should  be  used. 

The  maximum  base  displacement  for  the  optimized  parameters  was  also 
decreased,  but  the  large  oscillations  seen  in  the  first  example  did  not  occur.  Figure  21 
shows  the  base  displacements  for  the  initial  and  optimized  designs. 

Figures  22  and  23  show  hysteresis  loops  for  the  initial  and  optimized  designs, 
respectively.  Once  again,  the  optimized  hysteresis  loop  has  much  more  area  enclosed 
than  the  initial  design,  indicating  much  more  of  the  earthquake  energy  was  absorbed,  and 
resulting  in  the  decreased  top  floor  acceleration  magnitude  as  well  as  the  decreased  base 
displacement. 
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Isolator  Displacement  (in) 


Figure  22.  Hysteresis  Loop  for  Initial  Bilinear  Model  Parameters 


Figure  23.  Hysteresis  Loop  for  Optimized  Bilinear  Model  Parameters 
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VIII.  CONCLUSION  AND  RECOMMENDATIONS 


A.  CONCLUSION 

Use  of  the  RBBIF  algorithm  in  an  iterative  optimization  scheme  is  very 
appealing.  The  RBBIF  algorithm  is  very  fast,  computationally  efficient,  and  capable  of 
determining  the  transient  response  of  multi-story,  nonlinear,  base  isolated  structures.  The 
objective  of  this  thesis  was  to  utilize  this  algorithm  in  developing  an  optimization  scheme 
which  detennines  parameters  for  a  library  of  linear  and  nonlinear  seismic  isolator  models. 
By  also  utilizing  the  MATLAB  Optimization  Toolbox,  the  computational  efficiency  of 
the  resulting  optimization  program  was  further  increased. 

The  program  is  also  very  versatile.  Numerous  items  of  design  interest  can  be 
incorporated  into  the  objective  and  constraint  functions.  Both  horizontal  and  vertical 
responses  can  be  obtained,  although  one-dimensional  analysis  is  usually  the  primary 
concern.  In  addition,  any  design  earthquake  can  be  used  in  the  program. 

The  examples  show  that  the  optimization  program  can  be  very  fast.  The  objective 
functions  in  the  examples  were  to  minimize  peak  magnitudes  of  base  displacement  and 
top  floor  acceleration  to  some  fraction  of  the  initial  maximum  response.  The  program 
accomplished  these  objectives  in  a  very  short  time.  However,  even  though  the  objective 
function  was  successfully  minimized  in  each  example,  parameters  were  different  for  each 
of  the  first  three  examples  based  on  the  initial  starting  point.  The  fourth  example  also 
showed  that  the  optimal  parameters  determined  from  the  user-defined  constraints  may  not 
actually  support  a  good  design. 


B.  RECOMMENDATIONS  FOR  FUTURE  WORK 

The  optimization  program  does  not  include  three-dimensional  motion,  so  further 
study  in  this  area  would  be  useful  in  order  to  validate  the  program  against  a  real  structure 
and  a  FE  model  of  that  structure.  Additionally,  the  uplift  isolator  is  currently  modeled  as 
a  linear  spring  per  [Ref.  4],  but  nonlinearities  can  be  incorporated.  There  have  been 
studies  on  two  and  three-dimensional  motion  to  determine  rocking,  rolling,  and  uplift 
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effects  on  base  isolators  [Ref.  22].  Infonnation  such  as  this  may  be  useful  in  developing 
a  nonlinear  uplift  isolator  which  can  be  used  for  further  validation  of  the  program. 
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